Fast adaptive elliptical filtering using box splines 
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Abstract 

We demonstrate that it is possible to filter an image with an elliptic 
window of varying size, elongation and orientation with a fixed computa- 
tional cost per pixel. Our method involves the application of a suitable 
global pre-integrator followed by a pointwise-adaptive localization mesh. 
HH We present the basic theory for the ID case using a B-spline formalism 

^3 and then appropriately extend it to 2D using radially-uniform box splines. 

O The size and ellipticity of these radially-uniform box splines is adaptively 

controlled. Moreover, they converge to Gaussians as the order increases. 
Finally, we present a fast and practical directional filtering algorithm that 
has the capability of adapting to the local image features. 

oo 1 Introduction 

CO 

qq" The most common smoothing operator is the Gaussian filter. For that reason, 

f^*) it is of practical interest to design efficient directional filtering strategies based 

ON on this type of filter. Fast recursive solutions for space-invariant Gaussian-like 

filtering have been developed [1] but the space- variant ones are subject of current 
^ research. To date, two algorithms have been proposed: one that works in ID and 

uses B-splines for fast computations of the Continuous Wavelet Transform [3], 
and the other proposed in Computer Graphics that essentially does a rectangular 
smoothing using repeated integration [2]. 

In this paper, we propose a general technique for iV-directional adaptive fil- 
tering using box spline formalism and discuss its implementation for the special 
four-directional case. Although the derivation is rather involved, the final so- 
lution is quite simple (convolution-like with an adaptive mesh) and rather easy 
to implement (cf. Eqn. (15)). Also, the algorithm has a constant computa- 
tional cost per pixel as the support of the adaptive mesh is independent of the 
pointwise-adaptive scale-vector. 

The paper is organized as follows. In § 2, we revisit B-splines and certain 
associated linear operators. In § 3, we describe an efficient B-spline based adap- 
tive filtering technique, initially proposed in [3]. We then introduce the family 
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of radially-uniform box splines and some related linear operators in § 4 and 
§ 5, respectively. Finally, in § 6, we propose the adaptive directional filtering 
strategy followed by the description of a fast algorithm in § 7. 

2 ID Linear Operators and B-splines 

We first introduce two linear, shift-invariant operators, initially defined for real- 
valued functions f(x) and then appropriately applied to real-valued sequences 1 

g[k}. 

Definition 2.1 The finite- difference (FD) operator A™ of order n e Z + and 
scale a € R+ is specified as 

n 

A n J(x) = J2d n a [k]f(x-ak) (1) 
where d%[k] = a-™(-l) fe ("), forO<k<n and 0, else. 

In the Fourier domain, we have: A™/(w) = A™(e ja0J )f(w), with A"(e ja ") = 
J2k d™[k]e~ ]aujk = a~™(l — e -J aw )" interpreted as the (27r/a)-periodic frequency 
response of the FD filter d™ [k] . 

When acting on sequences g [k] , the FD filter acts as a discrete convolution 
operator. For integer a, we have A™g[m] = J2k=o da[k]g[m — ak]; for non-integer 
a, some form of interpolation is necessary as g is not defined for non-integer 
arguments. 

Definition 2.2 The running-sum (RS) operator A^ 1 of scale b G R+ is given 
by 

OO 

A-\f(x) = bJ2.f(x-bk). (2) 

fc=0 

For integer b, the RS operator applied to sequences g[k] corresponds to the 
digital filter u\j[k] = b, for k € M\o and 0, otherwise, with the correspondence: 
Afc~ 1 5[ n ] = Sfcez ^bl^M 71 — fc]- Note that y — Ub * g can be implemented very 
efficiently using the recursive equation y[m] = y[m — b] + bg[m], with appropriate 
boundary conditions [3] . The n-fold composition of the above RS operator will 
be denoted by A~" with frequency response A^ n {e jbu} ) = b n (l - e'^)- 71 ; 
for integer b, the corresponding filter u^[k] is specified by ^2 keZ [k]e ~i° jk = 
b n (l - e-^Y n . 

Finally, we introduce the family of symmetric B-splines [4] that are closely 
related to the above linear operators; it is the following Fourier domain definition 
that makes the link apparent. 

J By sequences we will mean functions denned on the Cartesian lattice Z d , where d is the 
dimensionality of the signal. 
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Definition 2.3 The symmetric B-spline f3™(x) of degree n G No and of scale 
a G R+ is specified by the Fourier transform 

1 / P ja.u/2_ -jaw/2 

«C)-iP«( * ) < 3 > 

Specifically, B-splines of arbitrary scales can be expressed in terms of integer- 
scaled B-splines: 

(%(x) =[K +1 o^- 1 [n+1) )^{x + r) (4) 

with t = (a — l)(n + l)/2, using the FD and RS operators [3]. It is also worth 
mentioning that /3„(x) = 1/a, for x € (—a/2, a/2] and else; it is the piecewise- 
constant function rect(x/a). In the sequel, we simply denote it by (3 a {x). 



3 Scale- Adaptive Filtering 

In this section, we revisit our ID space-variant filter [3] based on the projections 
s[m] = (f(x),P^ 2 (x — m)),m G Z, of a continuous signal model f(x) of the 
discrete signal with scaled B-splines /3™ 2 . The scale a controls the degree of 
smoothing applied around each sample. 

Specifically, given the signal samples f[k], we consider the following B-splinc 
model f(x) = Xlfegz c [fc]/3 ni ( x — k) for the continuum, with the interpolation 
constraint f(x)\ x= k = f[k], fceZ. The expansion coefficients are then obtained 
by the digital filtering c = f * (6™ 1 ) -1 , where is the convolution inverse 

of the B-spline interpolation filter b ni [k] = (3 ni (k), k € Z [4]. It then turns out 
that the projections s[m] can be efficiently realized using (4): 

Proposition 3.1 The B-spline projection s[m] can be computed in two steps: 

1. N on- Adaptive Step: The running sum filter t4" 2+1 ^ is applied to all the 
filtered sequence c[k] to get the integrated sequence 

g[m] = ^4" 2+1) [fc]c[m- k]. 
kez 

Then, the continuous domain pre-integrated signal can be written as F(x) — 
Ekez9[k]P ni+n2+1 (x + T-k). 

2. Adaptive Step: At each sample position, m G Z, and corresponding to 
a specific scale a, the integrated sequence g is then filtered using the FIR 
localization mask w[k] = A , Q l2+1 /3" 1+ ™ 2+1 ( fc + T ), k ^ z to 9 et the local 
projection s[m] — A"2 2+1 F(m) = X^fcez gW\w[m — k]. 
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4 Radially-Uniform Box Splines 



We now extend these ideas to 2D, where the additional feature of directionality 
needs to be addressed. Particularly, we devise the following tensor product 
between a ID smoothing operator and the identity convolution operator S : 
i\) i ► V'(O), both operating along orthogonal directions. 

Definition 4.1 The generating kernel ip a of scale a £ R + is defined as 

f a (x) = (3 a (x 1 )S(x 2 ), x = (x 1 ,x 2 ) e R 2 (5) 
The rotated versions ip a ,e(x) of the generating kernel are then obtained as 
<p a ,e( x ) = R 9 tpa(x) = (3 a (rJx)S(r'^x) 

where R g is the rotation operator defined as Rgf(x) — f(Rjx) via the rotation 
matrix 

„ / cos 9 — sin 9 \ / , \ 
y sin 9 cos 9 J \ ° J 

The most elementary box spline, the rect(x) function, can be generated by 
convolving two quadrature tensor B-splines; specifically, rect(x) = <pi i o(x) * 
fi,Tr/2{x) as per the above formalism. We generalize this idea to construct a 
family of box splines, which we call the radially-uniform box splines. 

Definition 4.2 The radially-uniform box spline is defined to be 

Pa(x) = ( i Pa 1 ,e 1 *---*<Pa N ,e N )(x) (6) 

where N £ N,iV > 2 is the directional order, a — (ai, . . . , ajv) € R+ is the 
scale vector, and 6k = (k — 1)tt/N, 1 < k < N are the rotation angles. 

Intuitively this means that we can construct 2D box splines by convolving arbi- 
trary number of tensor B-splines, rescaled by arbitrary amounts but uniformly 
distributed radially. In retrospect, rect(x) = 0k ^(x). 

Importantly, the covariance (moment) of these smoothing kernels can be ar- 
bitrarily controlled by suitably choosing the order TV and the pairs {{a-k,8k)}k=i- 
The remarkable fact is that as the order increases, these box splines become more 
Gaussian-like. In particular, we have the following result: 

Theorem 4.3 Let {Pa(N)} N > 2 be a sequence of box splines corresponding to 
the the scale-vector sequence {a(iV)}jv>2 with components given by ak(N) = 
CTy / 24/]V, 1 < k < N. Then we have the following convergence 

^«<»)W=2^ eX "(- ! Sf)- < 7 » 

Yet another form of convergence is achievable based on the iterated convolu- 
tion of the box spline (3^ , of a specific order N, with itself. Based on the Central 
Limit Theorem, existence of a sequence of iterated box splines that converges 
to a Gaussian can also be demonstrated. 
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5 2D Linear Operators 



We now extend the definitions of the FD (1) and the RS (2) operator to 2D; 
the key aspect is to preserve (4) in relation to the radially-uniform box-splines. 

Definition 5.1 The FD operator A a ^g of scale a g R + operating in the direc- 
tion 9 g [0, 7r) is specified by 



A a>e f(x) = ^f(x)-f(x-arg)). 



The TV-directional FD operator A^ corresponding to the scale-vector a = 
(ai, . . . , ajv) is then defined as the composition A^ = A ai ,0i o ■ ■ ■ o A aNt g N , 
where {A Uk ,e k }^ =1 are FD operators of scale Ofe g R+ operating in the direction 
Ok = (k — l)ir/N, repectively. 

Definition 5.2 The RS operator A^l of scale b g R+ operating in the direction 
9 g [0, 7r) is given by 

oo 

^,lf{x)=b^f{x-kbr ). (8) 

fe=0 

The TV-directional RS operator A^ N corresponding to the scale- vector b = 
(6i, . . . , b N ) is then defined as A b N = A^o- • ■oAj^, where {A^ e Jf =1 are 
RS operators of scale bk operating in the direction 9k = (k—l)ir/N, respectively. 

Importantly, based on (8), the application of A^ on a discrete sequence 
g[k] is then (non-uniquely) given by 

oo oo N 

A b N g[m] = &i ^ ■ ■ ■ 6jv ^ fint ( m -J2 k 3 b i r 0i) ( 9 ) 

fe 1= fe N =0 j=l 

where gi n t(x) is some form of interpolation of the discrete sequence g[k]. The 
good news is that for a specific choice of TV and b, the operator A^ N admits 
a convolution kernel (filter) and (9) then has a unique form; this motivates the 
following definition and subsequent simplification. 

Let < 9 < 7r and b g R+ be such that 6cos#, &sin6» g Z. Then the RS 
filter Ub,e[fc] of scale b and acting in the direction 9 is defined as 

. 1 6, for fc g (bcos9, 6sin#)N 
^ [fe] = \0, else (10) 

The TV-directional RS filter [fc] corresponding to the integration scale- vector 
b = . . . is then defined as the 2D convolution u b = u^fa * ■■■ * 
u b N ,e N of the RS filters {ub k ,e k }k=i °f sca l e b k and acting in the direction 9 k ; 
provided all the constituent filters are well-defined. Then (9) is uniquely given 
by A b N 9[m] = T, keZ 2 u£ [k]g[m - fc]. 
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Importantly, as in the ID case, it turns out that (3% can be factorized using 
the above linear operators as follows 

^(^(AfoA^^ + r) (11) 

where r = 0.5^^^ =1 (a^ — &&) cos 6^, X)tLi( a fc — &fc) sm ^fc) is the shift-vector. 

The role of the scale-vector b in (11) is to tie the RS operator A^ N to the 
underlying lattice Z 2 , and will be used to advantage in the sequel. 

6 Space- Variant Adaptive Filtering 

Our goal is to adaptively filter the signal f(x), obtained by interpolating the 
2D signal samples /[fe], with appropriately elongated and orientated box splines 
0a (x) by suitably adjusting the scale- vector a at the different spatial locations. 
In other words, given some pre-assigned scale- vector map a : Z 2 — > R+, we 
want to compute the projections s[m] = (f,Pa(m)(' ~ m ))' at spatial locations 

m e z 2 . 

First, as in the ID case, we consider the continuous representation f{x) = 
J2c[k] <p{x — fe), where <j>{x) is a 2D interpolating function. The expansion 
coefficients are given by c — f*{b n )~, where (6") _1 is the convolution inverse of 
the interpolation filter b n [k] = </>(k),k G Z 2 . Using (11), the desired projections 
can then be efficiently realized in two stages: 

(1) Non-Adaptive Step: The entire image / is filtered once using the fixed- 
scale RS filter to give the pre-integrated image: 



,A h N c\m], with interpolation 
9b[m} = { b N 1 \ ...... , .. 12 

{u b * c) [m] , without interpolation 



This also gives us the continuous pre-integrated image F(x) = J2kez 2 9b[k](4> 
^){x + T-k); 
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(2) Adaptive Step: The pre-integrated image gb is then locally filtered us- 
ing the pointwise localization operator &a(m) to §^ ve tn e adaptively filtered 
output 

s[m] = A^ m) F(k) = J2 ^a (m) [k}g b [m - k] (13) 
fcez 2 

at each location m G Z 2 corresponding to the scale-vector a(m), where w ( m ) [fe] = 
A^ Tfi ^*/3^ r (fc+T), feG Z 2 , is the localization mask, and r = 0. 5(^^ =1 (afc(m) — 
bk)cos9k,J2k=i( a k( m ) ~ bk)sin9k) is the pointwise shift-vector. 

7 Fast Elliptical Filtering Algorithm 

Finally, we focus on the family of four-directional box spline j3^(x) — (<p ai ,o * 
^Pa 2 ,-K/A * ( Pa 3 ,7r/2 * Vaz ,3-ir/i) ( x ) > corresponding to the directional order N = 4 
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Figure 1: Affine Mesh Geometry: The pair (•,•) denotes the position 
of a mesh vertex and the corresponding weight. 

and scale-vector a — (a%, . . . ,(14) in (6). It is interesting to note that 2(3"^(x), 
corresponding to a — (1,^/2,1,^/2), is also known as the ZP (Zwart-Powell) 
element [5] in box spline literature. 

Below, we outline the corresponding implementation aspects: 
(1) Non-Adaptive Step: First, we simplify the pre-filtcring by selecting the 
kernel <fi(x) = S(x), with the result that c[k] = f[k]. Further, the fact that 
(12) can be implemented (without interpolation) using the RS filter U(*[fe], with 
b = (1, y/2, 1, y/2), further simplifies the computation. In particular, the pre- 
integrated image can be expressed as 

g b [k\ = (ui, * "72,^/4 * M l,^/2 * u y/2,3n/i * ^ ^ ( 14 ) 

Now, due to this tensor structure, (14) can be then efficiently implemented in 
a recursive fashion in four steps, namely 

(1) Computation of Fq = u%fi * / using Fo[ki, fa] = f[fa, fa] + Fo[k% — 1, fa]. 

(2) Computation of = * Fq using F w /±[fa,fa] = V2F [fa , fa] + 

F„/4,[fa -l,k 2 - 1]. 

(3) Computation of F^/2 = "i.tt/2 * F^/ 4 using F w / 2 [fa,fa] = F w /±[fa, fa] + 

F„ /2 [fa,fa-l]. 

(4) Computation of g b = 37T / 4 *F n/2 as 9b[fa,fa] = V2F n/2 [fa,fa]+gb[fa + 

hfa-i}; 
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(2) Adaptive Step: Finally, simplifying (13), we get 

16 

s [ m ] = ^W]9i( m + T - x i) (15) 

i=i 

where gi(x) = gb[k][3^(x — k), is the ZP interpolation [5] of the discrete 
sequence gb[k], and f)[«] = (— l) 4+1 a, 1 < i < 16, are the non-zero weights of the 
affinc FD mesh, with a = (a\a 2 asa4)~ 1 . 

The corresponding (relative) positions {xi}j=\ C R 2 of the mesh vertices 
are shown in Table 1, with the convention a'j = aj/\/2, for j = 2,4. Figure 
1 gives the spatial ordering of the mesh vertices. The shift-vector in (15) is 
specified as r = (t\,t 2 ), where T\ — (\f2a\ + a 2 — a 4 — \/2)/2^/2 and r 2 = 
(a 2 + V2a 3 + C14 — 3 v / 2)/2v / 2- Note that r, f)[z] and £Ci are in fact defined 
pointwise in (15) using the scale- vector map a : Z 2 — > (cf. (13)); we 
dropped the pointwise index m just to simplify the equation. 



Table 1: Mesh Vertices 



«i : (0,0) 


#9 : (a x + a' 2 - a' 4 , a 3 + a' 2 + a' 4 ) 


x 2 : (oi,0) 


x w : (a 2 - a' 4 , a 3 + a' 2 + a' 4 ) 


^3 : (01,03) 


Xn ■ (a'2 - 04, a 2 + 04) 


aJ4 : (0,a 3 ) 


X12 ■ (ai + a' 2 - a 4 , a 2 + a 4 ) 


a; 5 : (a x + a' 2 , a' 2 ) 


®13 : (-a 4 ,a 3 + a 4 ) 


a; 6 : (01 + a' 2 , a 3 + a' 2 ) 


x 14 : (-a' 4 ,a' 4 ) 


x 7 : (a' 2 ,a 3 + a' 2 ) 


x 15 : (ai - a 4 ,o 4 ) 


x 8 : (02,02) 


a3i6 : (01 - a' 4 , a 3 + a 4 ) 



Note that the algorithm has a fixed computational cost per output pixel as 
the size of the support of the FD mask in (15) is independent of the scale- vector. 
Specifically, the number of non-null weights is 4 x 4 = 16, i.e., 4 clusters of 4 
points each, as shown in Figure 1. 

8 Conclusion 

We presented a novel elliptical filtering algorithm with fixed computational cost 
per output pixel. Our main goal was to formalize the adaptive elliptical filtering 
strategy using box splines and propose a fast and practical implementation 
algorithm. Computation of the scale-vector map presents a separate challenge 
in itself and will be discussed elsewhere. 
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